{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multipolar expansion approach for solving dipole radiation problems in presence of a dielectric waveguide\n",
    "\n",
    "Following Jackson's *lassical Electrodynamics*  book (chap 16, 1975 edition) as well as Prasad & Guo's 1997 paper (Optics Communications 136(1997) 447-460), I am posting some of the main result here regarding the classical radiation problem with a square waveguide. \n",
    "I am using Green's function formalism to formulate the radiated electric field from a dipole. \n",
    "The total field is the sum of a series of multiple scattering fields.\n",
    "Normally, one difficulty to calculate the scattered field is to avoid the sigularities in the integrations.\n",
    "That is why we are considering to use the multipolar expansion technique to solve the problem.\n",
    "\n",
    "As a summary, the Green's function tensor in free space responding at $\\mathbf{r}$ from a dipole source at $\\mathbf{r}'$ is given by\n",
    "$$\\mathbf{G}_0(\\mathbf{r}-\\mathbf{r}')=\\frac{1}{4\\pi}\\left(\\mathbb{1}+\\frac{1}{k_0^2}\\nabla\\nabla \\right)\\frac{e^{ik_0|\\mathbf{r}-\\mathbf{r}'|}}{|\\mathbf{r}-\\mathbf{r}'|}.$$\n",
    "The total electric field at position $\\mathbf{r}$ can then be given by\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathbf{E}(\\mathbf{r})&=\\mathbf{E}_{\\mathrm{inc}}(\\mathbf{r})+ \\sum_{n=1}^\\infty \\left[k_0(\\eta^2-1) \\right]^n\\int d^3r_1\\int d^3r_2\\cdots\\int d^3r_n \\\\\n",
    "&\\quad \\quad \\quad \\mathbf{G}_0(\\mathbf{r}-\\mathbf{r}_1)\\cdot \\mathbf{G}_0(\\mathbf{r}_1-\\mathbf{r}_2)\\ldots\\cdot \\mathbf{G}_0(\\mathbf{r}_{n-1}-\\mathbf{r}_n)\\cdot \\mathbf{E}_{\\mathrm{inc}}(\\mathbf{r}_n)\\\\\n",
    "&=\\mathbf{E}_{\\mathrm{inc}}(\\mathbf{r})+ \\sum_{n=1}^\\infty\\mathbf{E}_{\\mathrm{sc}}^{(n)}(\\mathbf{r}),\n",
    "\\end{align}\n",
    "$$\n",
    "where $\\eta^2=\\varepsilon/\\varepsilon_0$ and hence $\\eta^2-1=(\\varepsilon-\\varepsilon_0)/\\varepsilon_0=\\delta\\varepsilon/\\varepsilon_0$.\n",
    "\n",
    "In these notes, we are particularly interested in the local Green's function tensor responded at the dipole position in order to calculate the modified dipole decay rates in presence of a square waveguide. Therefore, we rewrite the $n$-th order scattering field at $\\mathbf{r}$ as\n",
    "$$\\begin{align}\n",
    "\\mathbf{E}_{sc}^{(n)} &= k_0^{2n}\\int dr_n\\mathbf{G}_n(\\mathbf{r},\\mathbf{r}_n)\\cdot \\mathbf{E}_{\\mathrm{inc}}(\\mathbf{r}_n)\n",
    "\\end{align}$$\n",
    "and define the $n$-th order scattering Green's function tensor \n",
    "$$\\begin{align}\n",
    "\\mathbf{G}_n(\\mathbf{r},\\mathbf{r}') &= \\int d^3r_1\\int d^3r_2\\cdots\\int d^3r_{n-1} \n",
    "\\mathbf{G}_0(\\mathbf{r},\\mathbf{r}_1)\\cdot \\frac{\\delta\\varepsilon(\\mathbf{r}_1)}{\\varepsilon_0}\\cdot \\mathbf{G}_0(\\mathbf{r}_1,\\mathbf{r}_2)\\cdot \\frac{\\delta\\varepsilon(\\mathbf{r}_2)}{\\varepsilon_0}\\ldots\\cdot \\mathbf{G}_0(\\mathbf{r}_{n-1},\\mathbf{r}')\\cdot \\frac{\\delta\\varepsilon(\\mathbf{r}')}{\\varepsilon_0}\\\\\n",
    "&=\\left(\\eta^2-1 \\right)^n\\int d^3r_1\\int d^3r_2\\cdots\\int d^3r_{n-1} \n",
    "\\mathbf{G}_0(\\mathbf{r},\\mathbf{r}_1)\\cdot \\mathbf{G}_0(\\mathbf{r}_1,\\mathbf{r}_2)\\ldots\\cdot \\mathbf{G}_0(\\mathbf{r}_{n-1},\\mathbf{r}')\n",
    "\\end{align}$$\n",
    "with all the integrations over the waveguide volume.\n",
    "Of course, the dipole source is at $\\mathbf{r}'$ position, but the scattering events happen through $\\mathbf{r}_1$ through $\\mathbf{r}_{n-1}$ positions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we use the spherical harmonics expansion that\n",
    "$$\\frac{e^{ik_0|\\mathbf{r}_i-\\mathbf{r}_j|}}{|\\mathbf{r}_i-\\mathbf{r}_j|}=4\\pi ik_0\\sum_{l=0}^\\infty\\sum_{m=-l}^l h_l^{(1)}(k_0r_{ij}^>)j_l(k_0r_{ij}^<)Y_{lm}(\\Omega_i)Y^*_{lm}(\\Omega_j),$$\n",
    "where $r_{ij}^>$($r_{ij}^<$) is the larger (smaller) of $r_i$ and $r_j$. \n",
    "\n",
    "Below, we will only consider the case that $\\mathbf{r}=\\mathbf{r}'$. \n",
    "It is well-known that the imaginary part of the $0$-th order Green's function tensor responded at the source can be given by\n",
    "$$\\begin{align}\n",
    "\\mathrm{Im}\\left[\\mathbf{G}_0(\\mathbf{r},\\mathbf{r})\\right] &=-\\frac{2}{3}k_0\\mathbb{1}.\n",
    "\\end{align}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Assume the origin is on the waveguide axis, and $\\mathbf{r}$ is always outside of the waveguide region. The first-order terms can be calculated as follows.\n",
    "$$\\begin{align}\n",
    "\\mathbf{G}_1(\\mathbf{r},\\mathbf{r}) \n",
    "&=\\left(\\eta^2-1 \\right)\\int_{WG} d^3r_1 \n",
    "\\mathbf{G}_0(\\mathbf{r},\\mathbf{r}_1)\\cdot \\mathbf{G}_0(\\mathbf{r}_1-\\mathbf{r})\\\\\n",
    "&=\\frac{\\eta^2-1}{(4\\pi)^2}\\int_{WG} dr_1 (4\\pi ik_0)^2\\left[\\left(\\mathbb{1}+\\frac{1}{k_0^2}\\nabla\\nabla \\right) \\sum_{l=0}^\\infty\\sum_{m=-l}^l h_l^{(1)}(k_0r)j_l(k_0r_1)Y_{lm}(\\Omega)Y^*_{lm}(\\Omega_1)\\right] \\\\\n",
    "&\\quad \\quad \\left[\\left(\\mathbb{1}+\\frac{1}{k_0^2}\\nabla_1\\nabla_1 \\right)\n",
    "\\sum_{l'=0}^\\infty\\sum_{m'=-l'}^{l'} h_{l'}^{(1)}(k_0r)j_{l'}(k_0r_1)Y_{l'm'}(\\Omega_1)Y^*_{l'm'}(\\Omega)\\right]\\\\\n",
    "&= -(\\eta^2-1)k_0^2 \\left[\\left(\\mathbb{1}+\\frac{1}{k_0^2}\\nabla\\nabla \\right) \\int_{WG} dr_1  \\sum_{l=0}^\\infty\\sum_{m=-l}^l h_l^{(1)}(k_0r)j_l(k_0r_1)Y_{lm}(\\Omega)Y^*_{lm}(\\Omega_1)\\right]\\\\\n",
    "&\\quad -(\\eta^2-1) \\int_{WG} dr_1 \\left[\\left(\\mathbb{1}+\\frac{1}{k_0^2}\\nabla\\nabla \\right) \n",
    "\\sum_{l'=0}^\\infty\\sum_{m'=-l'}^{l'} h_{l'}^{(1)}(k_0r)j_{l'}(k_0r_1)Y_{l'm'}(\\Omega_1)Y^*_{l'm'}(\\Omega)\\right]\\\\\n",
    "&\\quad\\quad\\quad\\quad\\quad\\quad\\quad\\quad \\left[ \\nabla_1\\nabla_1 \\sum_{l'=0}^\\infty\\sum_{m'=-l'}^{l'} h_{l'}^{(1)}(k_0r)j_{l'}(k_0r_1)Y_{l'm'}(\\Omega_1)Y^*_{l'm'}(\\Omega)\\right]\n",
    "\\end{align}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Given the square/rectangular shape of boundary condition, a spherical mode expansion of the free space Green's function solution is not a good choice as can be seen from the integrations above.\n",
    "It might be helpful to consider to use Hermite-Gaussian or Legarre-Gaussian modes to expand the free-space solution in order to obtain a simple form of integrals.\n",
    "For this purpose, the following references might be helpful: G. Beylkin *et al*. Journal of Computational Physics 228 (2009) 2770–2791 and F. Vico *et al.* Journal of Computational Physics 323 (2016) 191–203.\n",
    "\n",
    "Alternatively, one can expand the free-space Green's function into the Fourier space by using the Weyl identity \n",
    "$$\\frac{e^{ik\\sqrt{x^2+y^2+z^2}}}{\\sqrt{x^2+y^2+z^2}} = \\frac{i}{2\\pi}\\int\\int_{-\\infty}^{\\infty}\\frac{e^{ik_xx+ik_yy+ik_z|z|}}{k_z}dk_xdk_y$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.5.1",
   "language": "julia",
   "name": "julia-0.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
